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Abstract: In order to improve the convergence of multi-objective differential evolution (DEMO) algorithm 
while ensuring well distribution, a new method of center mutation-based DEMO (CM-DEMO) is proposed. 
Firstly, the form of mutation is improved, the center of the population in the current generation is taken as a 
base vector, and then the direction of difference vector is determined according to the fitness value of the three 
random vectors of individuals, secondly, the strategy of adaptive crossover probability is given, the crossover 
probability is determined according to the distribution of fitness value in the population. Test of benchmark 
functions show that CM-DEMO algorithm has faster convergence rate. Finally, CM-DEMO is applied to 
environmental economic dispatch of power system. Compared with other methods, the simulation results 


obtained demonstrate the feasibility and effectiveness of the proposed algorithm for solving the problem. 


Keywords: Differential evolution; Multi-objective optimization; Center mutation; Environmental economic 


dispatch 


1 Introduction 


Many real world problems can be formulated as optimization problems with multiple objectives. Since the 
first attempt to solve multi-objective optimization problems by using evolutionary algorithms, multi-objective 
evolutionary algorithms(MOEAs)have been much researched and are widely used to solve numerous 
applications in recent years (Coello, 2006; Nam & Park, 2000). MOEAs benefit from the evolutionary 
algorithm’s ability to generate a set of solutions concurrently in a single run thereby yielding several trade-off 


solutions. 


Differential evolution (DE) is a new generation evolutionary algorithm (EA) and has been successfully 
applied to solve a wide range of optimization problems. Differential evolution (DE) is a new type of 
evolutionary algorithm proposed by Storn and Price in 1995 (1997). It is simple yet powerful, and has 
been successfully used in solving single objective optimization problems. In recent years, some researchers 
have extended it to deal with multi-objective optimization problems, such as Pareto differential evolution 
(PDE) algorithm (Abbass, 2002; Xue, Sanderson & Graves, 2003), Pareto-based multi-objective differential 
evolution (PMODE) (Madavan, 2002) differential evolution for multi-objective optimization (DEMO) and 
adaptive differential evolution algorithm (ADEA) (Robic & Filipic, 2005; Qian & Li, 2008). 


However, most of these DE based multi-objective optimization algorithms suffer from premature convergence 
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at different degrees (Madavan, 2002; Robic & Filipic, 2005; Qian & Li, 2008). In this paper, we present 
a new multi-objective optimization of center mutation-based DEMO (CM-DEMO). The CM-DEMO has 
two improvements: a mutation operator composed of the modified base vector and differential vectors, the 
former is set as the center of all target vectors, and the latter is determined by the function fitness value of 
three randomly selected vectors; an adaptive crossover probability process according to the distribution of the 
function fitness value. For illustration, CM-DEMO was applied to solve the environmental economic dispatch 
problem, which consists of four interconnected cascade hydro plants and a thermal plant. Simulated results 


demonstrate the feasibility and effectiveness of the proposed method. 


The rest of the paper is organized as follows: The DE principle which CM-DEMO is based on is briefly 
described in Section 2. Afterward, in Section 3, we present CM-DEMO for solving problem in details. 
Section 4 presents the application study of CM-DEMO to a practical environmental economic dispatch 


problem. Section 5 outlines the conclusions followed by acknowledgements. 


2 Overview of Differential Evolution 


DE algorithm is a population based algorithm using three operators; crossover, mutation and selection. 


Several optimization parameters must also be tuned (Storn & Pric, 1997). The approach uses a population 


P® that contains NP n-dimensional real-valued parameter vectors (named a ) in generation g. According to 
Storn and Price, DE’s strategy can be described as follows. 


2.1 Mutation 


For each target vector xf, i=1,2,...,NP a mutant vector ee is generated according to (Storn & Pric, 
1997): 


y£" =x, +(%, x) i=1,2,...,NP (1) 


where integers n, n and n are chosen randomly in the range [1, N, ], and are different from 
each other. The mutation parameter F ([0, 2]) is a real, constant, user-supplied parameter that controls the 


amplification of the differential vector (named V ). 
2.2 Crossover 


In order to increase the diversity of the perturbed parameter vectors, crossover is introduced. The target vector 


gtl 
i 


is mixed with the mutated vector, using the following scheme, to yield the trial vector u 
& Pric, 1997): 


, as follows (Storn 


ij? 


za ie if (random() < CR) or j = randomRange(1,n) (2) 


x;;, otherwise 


In the above, function random () generates a random number in [0, 1], CR (€ [0, 1]) is the crossover 
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gtl 


parameter, and the integer j is a randomly chosen index in {1, 2,..., n} that ensures candidate u? get at 


least one parameter from trial parameter vector v* * but not all from xf. 


2.3 Selection 


When selection operation is implemented, a knock-out competition is played between the target vector xf 
gt 


and its corresponding trail vector u?  . According to the value of fitness, the better one will be selected for the 


next generation. Assuming that the fitness value is to be minimized, the selection operation can be expressed 
as (Storn & Pric, 1997): 


i 


pee p if u£#™ better than xë (3) 


X*; otherwise 


3 An Enhanced Multi-objective Differential Evolution: CM-DEMO 


In this section, by analyzing the mutation and crossover operation in the process of algorithm, the center 
mutation-based DEMO (CM-DEMO) algorithm is proposed. First, the CM-DEMO algorithm will be 
described in details, and then several typical and widely used benchmark test problems are chosen to test 
CM-DEMO. 


3.1 Center mutation operator 


From (1) we can conclude that DE algorithm randomly selects two individuals to calculate the difference as 
the difference vector, here the direction of difference vector is ignored. In this condition, the search capability 
is improved to some extent, but the convergence rate of the algorithm is slow down. In this section an 
operator called the center mutation operator is proposed, that is, each individual by mutation is around the 


center of contemporary population. The specific form of mutation operator is described as follows 
g+ Ce *( yE g *( yE g (4) 
VE HCl +F*(X} -X )+F*(X}-X3) 


Where, Xo is the best individual in the three random individuals; Xı and X2 are the other two random 


individuals; C is the center of population. 


Formula (4) shows that for each mutation individual, it is calculated based on the fitness value of the three 
random individual, then starts from the group center toward the of best individual X.. Here the direction 
which X, point to is the direction of difference vector, making the improved mutation operator not only retain 
a certain degree of randomness, but also the guidance of certainty is considered. Thus, under the action the 


center mutation operator, the new algorithm will have a faster convergence speed. 
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3.2 A daptive crossover probability 


In the standard DE algorithm, fixed CR is commonly used, that is in the process of solving problems, the 
crossover probability of each individual is identified with a certain value. However, in solving practical 
problems, the value should be adjusted appropriately based on the iteration cycle and the function values of 
each individual, making the algorithm acted in line with the characteristics of problem. Here a method of 


adaptive crossover probability is presented in Equation 5 (take the minimize problem as the case) 


F-F F-F 
max | —— H, = |, F > F (5) 
F, I Pin Pins = Fh 
CR, = 
: iF, F-F 
min| =], wt , other 
Po =k min Pax ~ min 


Where CR; is the crossover probability of each individual; F;and F, are the function values of the individuals 


Xj and X:; Fmin and Fmax are the values of the best and worst individuals in current generation. 


Formula (5) shows that, when F;>Fio, the crossover probability should be increased, in order to generate more 
mutation individuals Vi in the new individuals, the adjustment strategy is: by comparing the proportion which 
F; to the current generation and the proportion which F, to Fi, the greater proportion is choose as the value 
of crossover probability of the target individual; when Fi<Fo, the crossover probability should be reduced, 
then individual X; take more components in the new individuals, the adjustment strategy is: by comparing 
the proportion which Fj; to the current generation and the proportion which F, to Fi, the smaller proportion 
is choose as the value of crossover probability. This method can effectively adjust the composition of the 
individual generated according to the function value of individuals in iterative process, and then the search 


performance of algorithm is improved. 


3.3 Outline of CM-DEMO 


According to the above description of the improved algorithm, the outline of CM-DEMO is described as 


follows: 


Table 1. The main procedure of e-ODEMO 


Algorithm: Main procedure of the proposed CM-DEMO 
Initialize population randomly; 
While the halting criterion is not satisfied do 
Calculate the fitness function; 
Mutation according to formula (4); 
Calculate CR according to formula (5), crossover according to formula (2); 
Choose according to formula (2); 


end while 
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3.4 Performance test 


The test problems for evaluating the performance of our methods are chosen based on significant past studies 
in multi-objective evolutionary algorithms. We chose four problems from benchmark to test CM-DEMO 
(Zitzler & Thiele, 2000; Deb, Thiele, Laumanns & Zitzler, 2005). 


For every test problem: a crossover probability CR was set to 0.5, scaling factor F was set to 0.5. In order 
to make the comparisons fair, the population size NP was set to 100 and the algorithm was run for 250 


generations. 


Figures 1-4 show the Pareto fronts obtained by CM-DEMO and the real Pareto fronts of four ZDT test 
problems. As can be seen, solutions obtained by CM-DEMO scale very well in terms of convergence and 


widely-distributed. 


+ MEM ~ real Pareto front * CMDEMD —real pareto front 
i 1 
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Figure 1. Pareto front for ZDT1 Figure 2. Pareto front for ZDT3 
* MEM real Pareto front + MO ——real Pareto front 
1 1 
0.8 0.8 
0.6 0.6 
È È 
0.4 0.4 
0.2 0.2 
0 0 
0 0.2 0.4 Fl 0.6 0. 8 0 0.2 0.4 Fl 0.6 0.8 
Figure 3. Pareto front for ZDT4 Figure 4. Pareto front for ZDT6 


Tables 2 present the mean (boldfaced font above) and variance (underside) of the values of the convergence 
and diversity metric averaged over 10 runs (Deb & Jain, 2002; Van Veldhuizen, 1999). Results of other 
algorithms are taken from the literatures((NSGA- II, SPEA2, PDEA, DEMO, ADEA and DEMO/parent). 
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Table 2. Statistics of results on convergence metric y and diversity metric A 


Metrics Problems NSGA-II SPEA2 DEMO 
(real-coded) 


Y ZDT1 0.033482 0.023285 0.001083 
0.004750 0 0.000113 

ZDT3 0.114500 0.018409 0.001178 

0.007940 0 0.000059 

ZDT4 0.513053 4.9271 0.001037 

0.118460 2.703 0.000134 


ZDT6 0.296564 0.232551 0.000629 

0.013135 0.004945 0.000044 

A ZDT1 0.390307 0.154723 0.325237 
0.001876 0.000874 0.030249 

ZDT3 0.738540 0.469100 0.309436 

0.019706 0.005265 0.018603 

ZDT4 0.702612 0.823900 0.359905 

0.064648 0.002883 0.037672 

ZDT6 0.668025 1.04422 0.442308 

0.009923 0.158106 0.021721 


PDEA 


m W W M -M 


/ 
0.298576 
0.000742 
0.623812 
0.000225 
0.840852 
0.035741 
0.473074 
0.021721 


ADEA 


0.002741 
0.000385 
0.002741 
0.000120 
0.100100 
0.446200 
0.000624 
0.000060 
0.382890 
0.001435 
0.525770 
0.043030 
0.436300 
0.110000 
0.361100 
0.036100 


PDEMO 


0.005800 
0 
0.021560 
0 
0.638950 
0.500200 
0.026230 
0.000861 
/ 


—_— s “~ “ “ ™. 


CM- 
DEMO 
0.000028 
0 
0.00017 
0 
0.002190 
0.000097 
0.000016 
0 
0.146883 
0.000062 
0.433462 
0.000453 
0.157800 
0.000305 
0.114263 
0.001769 


The results of convergence metric(see Table 2) for these three problems show that CM-DEMO achieves good 
convergence, which is rather better than PDEA, DEMO, and DEMO/parent, and much better than NSGA- II 
and SPEA2. On ZDT3, DEMO achieves comparable results to CM-DEMO, but on ZDT1 and ZDT6, CM- 
DEMO performances better. For another metric A , showed in Table 2 CM-DEMO achieves much better 


results on four test problems than the other algorithms referred here. 


ZDT4 is a difficult optimization problem with large number of local Pareto fronts that tend to mislead the 
optimization algorithm. In Table 2 we can see that NSGA- II, SPEA2, PDEA and MODE all have difficulties 
in converging to the true Pareto front. DEMO, CM-DEMO performs better than the other algorithms. 


Problem Tamaki is constrained test problem and DTLZ1 is high-dimension problem (M=3). Figures 5 and 6 


show the Pareto fronts obtained by CM-DEMO. 


Figure 5. Pareto front for Tamaki 
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Figure 6. Pareto front for DTLZ1 
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From Figures 5 and 6 we can see that CM-DEMO handled constraints and high-dimension problem well, also 
it converged to the true Pareto front accurately. After achieving good performance on test problems, next we 


will apply CM-DEMO to a practical problem of environmental economic dispatch of power systems. 


4 Case Sstudy 


In order to validate the proposed procedure, a test hydrothermal system is taken for case study (Naresh & 
Sharma, 1999). The system consists of a multi-chain cascade of four hydro plants and three thermal units. 


The details data of the system considered here are the same as in Ref. 
4.1 Hydrothermal scheduling problem 


The solution of environmental economic dispatch problem aims to minimize the operation costs of thermal 
power plant and contaminative gas emission simultaneously, while satisfying a series of equality and 
inequality constraints (Talaq, El-Hawary & EI-Hawary, 1994; Yalcinoz & Köksoy, 2007; Abido, 2003). 


4.1.1 Problem objectives 
(1) Economy objective 


The generator cost curves are represented by quadratic functions of real power generation by that unit, the 
total fuel cost can be formulated as: 
TN, TN, ; 
min f= > > LA (POI = > > [ay + by Pi te, (PY +d sin(e, (P - Pi) (6) 
i=1 


t=1 i=l 


=] i= 


where a,,,b,andc,, are the cost curve coefficients of the ith thermal unit, P,(¢) is the output power of the ith 


si? 


thermal unit at period £, d, e, are the valve-point effects coefficients of the ith thermal plant; P7™ is the 


i 


minimum output limit of the ith thermal plant. 
(2) Emission objective 


Since SO2 and NOx emissions are generally taken to be proportional to the generator’s fuel consumption, 
total emission function to represent SO2 and NOx emissions are used in this paper is the same form as that of 


the fuel cost function. Hence, the total emission of the hydro-thermal system can be expressed as: 


T N 
fl i 


min f,=min > > [a, +8, PO + ya P; (O +n, exp, PO] (7) 


| i=l 


where a, Bus Yus Nus and 6, are emission curve coefficients of the ith thermal unit. 
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4.1.2 System constraints 
The problem subjects to is subjected to the following equality and inequality constraints: 


(3) Power balance constraints 


N, Nn 
SPAS -B-P =0, t2 8 


where P,, is the load demand at period t, P, „the total transmission line losses at period t, P” 


i,t 


is the output 
power of the jth hydro plant at period t, and Nh is the number of hydroelectric plants. 


(4) Real power output limits 


<P), i=l., N, and Pi, $ Pi < Pho i=l,...N, (4) 


j,min j j,max ? 


1 „Max 


where, Ph ia and P' are the lower and upper generation limits of the jth hydroelectric plant, respectively. 


(5) Reservoir storage volumes limits 


V 


j.min 


<V,,<V, 


j,max ? 


j=1,.., N, (9) 


where yV, min and y, max are the minimum and maximum storage volume of the jth reservoir, respectively. 


(6) Hydro plant power limits 


Q; min < 0, < O sacs J = L.. N, (10) 


where Q j,min and O; max are the minimum and maximum water discharge rate of the jth hydroelectric plant, 
respectively. 


(7) Initial and terminal reservoir storage volumes 


Fa = yee, Vg = eae j = L,..., N, (11) 


bevi d PAPEN os 
where V;" and V;"" are the initial and final storage volume of reservoir j. 


(8) Water dynamic balance equation with travel time 
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N, 


> 12 
Vi = V + lf =. =S ga Cx F Ser) -At ( ) 


k=l 


where /,, is the nature inflow rate of the jth reservoir at period t, Nj is the number of upstream units 


directly above the jth hydroelectric plant, ee is the spillage of the jth reservoir at time ¢, and Ty is the water 
transport delay from reservoir k to /. 


(9) Hydroelectric generation relation 


Pre = Cy VY + Cay OpY + Cay Vig Qpe + Cay Vig Coy Qy + Ce, (13) 
$12 PEL 


where Cy, Cy, C3, Cy, Cy, and Cg are the power generation coefficients of jth hydroelectric plant, Q,, is the 


water discharge rate of jth reservoir at period t, V;, is the storage volume of jth reservoir at period t. 


4.2 Solution methodology based on CM-DEMO 


In this section, some method of proposed CM-DEMO for solving generation scheduling of this hydrothermal 
system is described in details. Especially, a suggestion will be given on how to handle constraints of the 


problem. 
4.2.1 Initialization 


In the initialization procedure, the population is initialized by creating NP solutions randomly. For all 
solutions in the initial population, each element is randomly generated within the feasible real power output 
range to the constraint. For the present problem, it is discharge rate of each hydro plant and the power 


generated by each thermal unit. Thus the population is initialized as follows: 


Oi. 505s ` SO "Qn, KOMETE On 7: 
S S S S S S 
Ri 1,237073 iro o Pao Pv EN 
The dimension y = [(N, +N) *T]. Py is randomly generated between Q, min and O max; and P, is 


randomly generated between P’ and Pi Generally, the newly generated individuals do not satisfy all 


j,min j,max 


the constraints and need to be modified by the constraints handling method, described next. 
4.2.2 Mutation, crossover and selection 


New values of water discharge rate and power generation are generated through mutation and crossover 
operation according to (8) and (2) respectively. Now, selection is performed by calculate the fitness values of 
the different individuals. The individuals in the current population are evaluated in the objective space and 
then assigned a scalar value known as fitness. Depending on the fitness values, individuals will be selected 
to form the new population. Individuals which have a low fitness value have the chance to be selected. It 


is worth mentioning that the constraint-handling approach implemented in this study is that the unfeasible 
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solutions are penalized by assigning a very high value for their fitness. 


4.2.3 Constraint handling 


Firstly, When population initialization, crossover has been implemented, the new generated solution may not 
satisfy equality constraints (1) and (5). At present, penalty method is the most popular constraints handling 
strategy for dealing with this equality constraint at present by using penalty function to punish the infeasible 
solution during the selection procedure to ensure the priority of feasible ones. However, this strategy may 


degrade the efficiency of the algorithm remarkably for it requires multiple runs to tune the penalty factors. 


Secondly, we focus on handling the output capacity limits (2) and water release limits (4) when the proposed 


CM-DEMO method is applied to solve environmental economic dispatch problem. 


Despite the popularity of penalty functions, they have several drawbacks among which the main one is that 
they require a careful fine tuning of the penalty factors that accurately estimates the degree of penalization 
to be applied as to approach efficiently the feasible region. In order to keep the advantages of the penalty 
function approach and overcome drawback of choice penalty factors, this paper apply an effective constraint 
handling method for DE, which does not require to set any additional parameters in comparison with 
the original DE. Therefore, in order to strike a balance between computational efficiency and constraints 
handling, the following selection strategy is adopted by the proposed CM-DEMO method to choose the better 


solution while considering the constraints violation during the selection operation: 


(1) If solution P1 is feasible and solution P2 is infeasible, then P1 is favored. 


(2) If both P1 and P2 are feasible, the Pareto-dominance based selection will be implemented to decide which 
one is better, the selection operation is modified as follows: 

1) If the candidate dominates the parent, replace the parent by the candidate. 

2) If the parent dominates the candidate, the candidate is discarded. 

3) If the candidate and parent are non-dominated with each other, a new population of the size between NP 


and 2NP is created, and then add the non-dominated candidates and parents into it. 


(3) If both P1 and P2 are infeasible, then the one with smaller constraints violation is favored. 


4.3 Simulation results 


With the parameter settings listed in Table 3, the proposed method coded by Microsoft Visual C++6.0 
language on a Pentium-4 2.0GHz-based processor computer is applied to solve the optimal generation 
scheduling of this hydrothermal system. The hourly each hydro plant power generation are showed in Table 4. 


The hourly each reservoir release and storage trajectories are showed in Figures 7 and 8 respectively. 
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Method Generation number 


Table 3. Parameter settings for CM-DEMO 


Population size 


CM-DEMO 1000 50 0.5 


Table 4. Hourly hydrothermal power generation scheduling (Unit 


Hour 1 2 3 4 5 6 7 8 9 10 

Pai 73.77 79.17 80.44 75.14 78.26 68.15 72.03 79.44 79.42 88.44 
Pn2 52.93 50.04 50.96 53.97 54.09 56.03 61.14 60.82 76.02 67.39 
Pas 34.60 24.82 27.47 3486 32.53 22.29 38.23 33.18 37.06 26.93 
Pha 131.8 | 129.0 125.7. 122.7 115.7 130.8 190.4 223.2 283.0 300.5 
Psi 104.9 162.6 174.9 104.4 107.9 174.9 174.7 103.9 174.9 175 

Pa 198.7 193.0 125.2 207.4 196.7 209.2 220.4 211.1 206.9 209.8 
Ps3 153.0 141.3 115.1 51.41 84.73 138.4 192.9 298.2 232.5 211.7 
Hour 13 14 15 16 17 18 19 20 21 22 

Pri 90.27 78.73 83.22 86.94 84.29 81.08 71.78 65.18 57.31 66.57 
Pra 82.86 61.42 70.39 72.20 72.34 75.35 63.85 69.08 82.69 88.51 
Phr 33.14 40.57 42.01 46.76 49.24 51.27 51.76 54.12 54.81 57.39 
Pra 288.6 304.0 282.5 302.2 295.0 293.2 295.4 276.4 295.3 259.3 
Por 174.4 174.9 174.2 174.9 174.9 175 162.2 | 111.1 | 137.9 | 123.2 
Ps 230.5 215.2 | 210.5 ; 228.9 210.8) 217.2 ; 1995 | 179.3 ; 126.7 | 125.9 
Ps3 210.0 155.0 147.0 147.9 163.2 | 226.8 225.2 | 294.6 155.1 139.0 


Crossover parameter Scaling factor 


0.1 


: 104KW) 

11 12 
84.86 84.38 
56.32 85.79 
37.47 28.91 
289.4 299.1 
174.8 174.8 
234.7 223.7 
222.2 253.1 
23 24 
59.97 59.23 
73.67 61.52 
56.71 57.91 
288.8 279.5 

125.72 | 120.75 
195.0 171.0 
50 50 


—eFeservoir 1 —=— Reservoir 2 —*—Reservoir 3 —<—Reservoir 4 


ay 
au oO 


Di schar ge( 10°nf/h) 
on 


© 


Figure 7. Hourly hydro plant discharge 
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Figure 8. Hourly hydro plant storage 


In the meantime, to validate the results obtained with the proposed CM-DEMO method, the same problem is 
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solved using DEMO, NSGA II, SPEA II methods, the results of test system are also summarized in Table 5 


for the convenience of comparisons. 


Table 5. Comparisons of schedule results with other methods 


Method Total fuel cost ($) Total emission (Ib) 
CM-DEMO 44781 16870 
DEMO 45050 17083 
NSGA II 45534 16974 
SPEA II 45431 16838 


It is clear from Table 3 that the total fuel cost and total emission obtained by CM-DEMO are much less 
compared to the corresponding values of the other methods. In addition, during these 20 independent 
simulations, it demonstrates that total fuel cost and total emission generate a variation in a small range with 
trial numbers when using the proposed CM-DEMO method. As can be seen from the simulation results of 
CM-DEMO method, the solutions are optimal and they also satisfy various constraints completely for solving 


environmental economic dispatch problem. 


5 Conclusions 


In this paper, a modified multi-objective differential evolution optimization algorithm based on center 
mutation mechanisms has been proposed. The algorithm replaces the original mutation and crossover 
operation to the improved form, which is center mutation and adaptive crossover. Simulation results and the 
comparison confirm the effectiveness and the superiority of the proposed approach over the other techniques 
in terms of the quality and precision of solution, so it provides an effective method to solve the environmental 


economic dispatch. 
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